﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace LPLib
{

    //public class PairwiseAligner
    //{
    //    /// <summary>
    //    /// The total score for the alignment as reported.
    //    /// </summary>
    //    public double AlignmentScore { get; protected set; }

    //    private double[, ,] scoreMatrix;
    //    private int[, ,] traceBackMatrix;

    //    private int dimQueryPhysical = -1;
    //    private int dimReferencePhysical = -1;
    //    private int dimQuery;
    //    private int dimReference;
    //}


    /// <summary>
    /// Provides the means to perform pairwise alignment.
    /// </summary>
    /// <typeparam name="T">The nucleotide type for the first polynucleotide.</typeparam>
    /// <typeparam name="U">The nucleotide type for the second polynucleotide.</typeparam>
    public class PairwiseAlignerAsymmetric<T, U>
    {
        /// <summary>
        /// The total score for the alignment as reported.
        /// </summary>
        public double AlignmentScore { get; protected set; }

        private double[, ,] scoreMatrix;
        private int[, ,] traceBackMatrix;

        private double[] cumulativeScoreQuery;

        private int dimQueryPhysical = -1;
        private int dimReferencePhysical = -1;
        private int dimQuery;
        private int dimReference;

        private AlignmentKernelAsymmetric<T, U> kernel;

        /// <summary>
        /// Sets the default scores for the first column.
        /// </summary>
        /// <param name="values">A vector of cumulative alignment scores taken as the alignment scores that occur if the polynucleotide local alignment begins on a row other than the zeroth.</param>
        public void SetFirstColumn(double[] values)
        {
            cumulativeScoreQuery = values;
        }

        /// <summary>
        /// Sets the default scores for the first column.
        /// </summary>
        /// <param name="value">The score to be accumulated over rows as the alignment scores that occur if the polynucleotide local alignment begins on a row other than the zeroth.</param>
        public void SetFirstColumn()
        {
            if (dimQuery < 1)
                return;

            cumulativeScoreQuery[0] = kernel.PriorQueryScore(0);
            for (int i = 1; i < dimQuery; i++)
            {
                cumulativeScoreQuery[i] = cumulativeScoreQuery[i - 1] + kernel.PriorQueryScore(i);
            }
        }

        /// <summary>
        /// Sets the default scores for the first column.
        /// </summary>
        /// <param name="value">The score to be accumulated over rows as the alignment scores that occur if the polynucleotide local alignment begins on a row other than the zeroth.</param>
        public void SetFirstColumn(double defaultValue)
        {
            if (dimQuery < 1)
                return;

            cumulativeScoreQuery[0] = defaultValue;
            for (int i = 1; i < dimQuery; i++)
            {
                cumulativeScoreQuery[i] = cumulativeScoreQuery[i - 1] + defaultValue;
            }
        }

        /// <summary>
        /// Constructor.
        /// </summary>
        /// <param name="_pair">The polynucleotide pair to be aligned.</param>
        /// <param name="_comparator">A comparator of the appropriate types for use in the alignment.</param>
        public PairwiseAlignerAsymmetric(AlignmentKernelAsymmetric<T, U> kernel)
        {
            this.kernel = kernel;
            if (kernel.GetQuery() == null) //  || _comparer.GetReference() == null)
            {
                dimQuery = dimQueryPhysical = 0;
            }
            else
            {
                dimQuery = dimQueryPhysical = kernel.GetQuery().Seq.Length;
            }

            if (kernel.GetReference() == null)
            {
                dimReference = dimReferencePhysical = 0;
            }
            else
            {
                dimReference = dimReferencePhysical = kernel.GetReference().Seq.Length;
            }

            // allocate matrices
            scoreMatrix = new double[dimQueryPhysical, dimReferencePhysical, 3];
            traceBackMatrix = new int[dimQueryPhysical, dimReferencePhysical, 3];
            cumulativeScoreQuery = new double[dimQueryPhysical];
            SetFirstColumn();
        }

        /// <summary>
        /// Constructor.
        /// </summary>
        /// <param name="kernel">A comparator of the appropriate types for use in the alignment.</param>
        /// <param name="_dim0">The initial first dimension of the score and traceback matrices.</param>
        /// <param name="_dim1">The initial second dimension of the score and traceback matrices.</param>
        public PairwiseAlignerAsymmetric(int _dim0, int _dim1, AlignmentKernelAsymmetric<T, U> kernel)
        {
            this.kernel = kernel;
            dimQuery = dimQueryPhysical = _dim0;
            dimReference = dimReferencePhysical = _dim1;

            // allocate matrices
            scoreMatrix = new double[dimQueryPhysical, dimReferencePhysical, 3];
            traceBackMatrix = new int[dimQueryPhysical, dimReferencePhysical, 3];
            cumulativeScoreQuery = new double[dimQueryPhysical];
       }

        public void SetQuery(Polymer<T> p, double[] defaultScores = null)
        {
            kernel.SetQuery(p);
            dimQuery = p.Seq.Length;
            if (dimQuery > dimQueryPhysical)
            {
                dimQueryPhysical = dimQuery;
                scoreMatrix = new double[dimQueryPhysical, dimReferencePhysical, 3];
                traceBackMatrix = new int[dimQueryPhysical, dimReferencePhysical, 3];
                cumulativeScoreQuery = new double[dimQueryPhysical];
            }

            if (defaultScores != null)
            {
                SetFirstColumn(defaultScores);
            }
            else
            {
                SetFirstColumn();
            }
        }

        public void SetReference(Polymer<U> p, double[] defaultScores = null)
        {
            kernel.SetReference(p);
            dimReference = p.Seq.Length;
            if (dimReference > dimReferencePhysical)
            {
                dimReferencePhysical = dimReference;
                scoreMatrix = new double[dimQueryPhysical, dimReferencePhysical, 3];
                traceBackMatrix = new int[dimQueryPhysical, dimReferencePhysical, 3];
            }
        }

        public double ScoreAlignedPair(AlignmentKernelAsymmetric<T, U> _comparer, PolymerPair<T, U> _pair)
        {
            // This method assumes that a query-reference pair has already been aligned and that a second reference, sufficiently similar to the first that a complete 
            // realignment is not necessary, is being substituted for the first.
            //if (cumulativeScore0 != null)
            //    SetFirstColumn();

            // for this method, we assume that _comparer.query = _pair.Polynucleotide1
            if (_pair.Index.Length == 0)
                return double.MinValue;

            int offset = Math.Min(_pair.Index[0][0], _pair.Index[0][1]);

            // firstx is the first position (in seqx coordinates) in the aligned _comparer.Pair where both sequences have nucleotides
            int first0 = _pair.Index[0][0] - offset;
            int first1 = _pair.Index[0][1] - offset;
            int[] lastIndexed = _pair.Index.Last();

            // the overhang is computed because the new allele may align to the query beyond the point where the original allele aligned.
            int overhang = Math.Min(_pair.Polymer0.Seq.Length - lastIndexed[0], _pair.Polymer1.Seq.Length - lastIndexed[1]) - 1;

            // analogously for lastx
            int last0 = lastIndexed[0] + overhang;
            int last1 = lastIndexed[1] + overhang;

            double unmatchedScore = 0;
            if (first0 > 0)
            {
                unmatchedScore = cumulativeScoreQuery[first0 - 1];
            }

            int newLength = _pair.Index.Length + offset + overhang;

            double[] newScore = new double[newLength];
            int[] traceback = new int[newLength];
            int[][] newIndex = new int[newLength][];

            newScore[0] = unmatchedScore + kernel.ComparisonScore(first0, first1);
            traceback[0] = 1;
            newIndex[0] = new int[] { first0, first1 };
            // compute scores up to and including the start of the existing index  
            for (int kNew = 1; kNew < offset + 1; kNew++)
            {
                int i0 = first0 + kNew;
                int i1 = first1 + kNew;
                newScore[kNew] = kernel.ComparisonScore(i0, i1);
                if (newScore[kNew - 1] > cumulativeScoreQuery[i0 - 1])
                {
                    newScore[kNew] += newScore[kNew - 1];
                    traceback[kNew] = 0;
                }
                else
                {
                    newScore[kNew] += cumulativeScoreQuery[i0 - 1];
                    traceback[kNew] = 1;
                }
                newIndex[kNew] = new int[] { i0, i1 };
            }

            // compute scores through the end of the existing index
            int endK = _pair.Index.Length + Math.Min(0, overhang);
            for (int k = 1; k < endK; k++)
            {
                int i1 = _pair.Index[k][1];
                if (i1 == _pair.Polymer1.Seq.Length) // if there is a gap in the original allele between the end of the original index 
                {                                   // and the end of the new allele, this condition will get triggered
                    newLength -= (endK - k);
                    break;
                }

                int kNew = offset + k;
                int i0 = _pair.Index[k][0];
                if (i0 > -1 && i1 > -1) // both sequences have a nucleotide at this position
                {
                    newScore[kNew] = kernel.ComparisonScore(i0, i1);
                    if (newScore[kNew - 1] > cumulativeScoreQuery[i0 - 1])
                    {
                        newScore[kNew] += newScore[kNew - 1];
                        traceback[kNew] = 0;
                    }
                    else
                    {
                        newScore[kNew] += cumulativeScoreQuery[i0 - 1];
                        traceback[kNew] = 1;
                    }
                }
                else if (_pair.Index[k][0] < 0) // sequence0 has a gap at this position
                {
                    // if this gap is aligned opposite a gap in sequence1, no penalty
                    if (kernel.IsGap(_pair.Polymer1.Seq[i1]))
                    {
                        newScore[kNew] = 0;
                    }

                    // note that the first index value cannot be negative for either sequence
                    if (_pair.Index[k - 1][0] < 0)
                    {
                        newScore[kNew] = kernel.ContinueDeletionScore(i0, i1);
                    }
                    else
                    {
                        newScore[kNew] = kernel.OpenDeletionScore(i0, i1);
                    }
                    newScore[kNew] += newScore[kNew - 1];
                    traceback[kNew] = 0;
                }
                else // sequence1 has a gap at this position.
                {
                    if (_pair.Index[k - 1][1] < 0)
                    {
                        newScore[kNew] = kernel.ContinueInsertionScore(i0, i1);
                    }
                    else
                    {
                        newScore[kNew] = kernel.OpenInsertionScore(i0, i1);
                    }
                    newScore[kNew] += newScore[kNew - 1];
                    traceback[kNew] = 0;
                }
                newIndex[kNew] = _pair.Index[k];
            }

            // continue past the end of the old index
            for (int deltaKNew = 0; deltaKNew < overhang; deltaKNew++)
            {
                int kNew = _pair.Index.Length + offset + deltaKNew;
                int i0 = lastIndexed[0] + deltaKNew + 1;
                int i1 = lastIndexed[1] + deltaKNew + 1;
                newScore[kNew] = newScore[kNew - 1] + kernel.ComparisonScore(i0, i1);
                traceback[kNew] = 0;
                newIndex[kNew] = new int[] { i0, i1 };
            }

            // do traceback
            int stop = newLength - 1;
            double cumScoreLast = cumulativeScoreQuery[_pair.Polymer0.Seq.Length - 1];
            double totalScore = newScore[stop] + cumScoreLast - cumulativeScoreQuery[last0];
            double maxTotalScore = totalScore;
            // we have to go back through the index here
            for (int kNew = newLength - 1; kNew > 0; kNew--)
            {
                int i0 = newIndex[kNew][0];
                if (i0 < 0)
                    continue;

                totalScore = cumScoreLast - cumulativeScoreQuery[i0] + newScore[kNew];
                if (totalScore > maxTotalScore)
                {
                    stop = kNew;
                    maxTotalScore = totalScore;
                }
            }
            int start = stop;
            while (traceback[start] == 0)
            {
                start--;
            }

            // form new _comparer.Pair index and new _comparer.Pair score array
            _pair.Index = new int[stop - start + 1][];
            _pair.Score = new double[stop - start + 1];
            Array.Copy(newIndex, start, _pair.Index, 0, stop - start + 1);
            Array.Copy(newScore, start, _pair.Score, 0, stop - start + 1);
            _pair.AlignmentScore = maxTotalScore;
            return maxTotalScore;
        }

        /// <summary>
        /// Fills the score matrix for pairwise alignment between the two Polynucleotides in the pair.
        /// </summary>
        public int[] FillScoreMatrix(bool qFivePrimeEndFree = true, bool rFivePrimeEndFree = true, bool qThreePrimeEndFree=true, bool rThreePrimeEndFree = true)
        {

            dimQuery = kernel.GetQuery().Seq.Length;
            dimReference = kernel.GetReference().Seq.Length;

            int i0 = 0;
            int i1 = 0;

            double scoreDiag, scoreDiagMax;
            double scoreDelete, scoreDeleteMax;
            double scoreInsert, scoreInsertMax;
            double jumpScore;
            bool bothFivePrimeEndsFree = qFivePrimeEndFree && rFivePrimeEndFree;

            // compute element [0,0]
            scoreMatrix[0, 0, 0] = double.MinValue; // cannot end on an internal gap
            scoreMatrix[0, 0, 1] = kernel.ComparisonScore(0, 0);
            traceBackMatrix[0, 0, 1] = -1;
            scoreMatrix[0, 0, 2] = double.MinValue; // cannot end on an internal gap

            // compute zeroth row
            for (i1 = 1; i1 < dimReference; i1++)
            {
                scoreMatrix[0, i1, 0] = double.MinValue;
                scoreMatrix[0, i1, 1] = kernel.ComparisonScore(0, i1);
                traceBackMatrix[0, i1, 1] = rFivePrimeEndFree ? -1 : 1;
                scoreMatrix[0, i1, 2] = double.MinValue;
            }

            // compute zeroth column
            for (i0 = 1; i0 < dimQuery; i0++)
            {
                scoreMatrix[i0, 0, 0] = double.MinValue;
                scoreMatrix[i0, 0, 1] = cumulativeScoreQuery[i0 - 1] + kernel.ComparisonScore(i0, 0);
                traceBackMatrix[i0, 0, 1] = qFivePrimeEndFree ? -1 : 1;
                scoreMatrix[i0, 0, 2] = double.MinValue;
            }

            // compute remaining rows
            for (i0 = 1; i0 < dimQuery; i0++)
            {
                for (i1 = 1; i1 < dimReference; i1++)
                {
                    // layer 0: insertions into sequence 0
                    scoreInsertMax = double.MinValue;
                    for (int iLayer = 0; iLayer < 2; iLayer++)
                    {
                        scoreInsert = scoreMatrix[i0 - 1, i1, iLayer] + (iLayer == 0 ? kernel.ContinueInsertionScore(i0, i1) : kernel.OpenInsertionScore(i0, i1));
                        if (scoreInsert > scoreInsertMax)
                        {
                            scoreInsertMax = scoreInsert;
                            traceBackMatrix[i0, i1, 0] = iLayer;
                        }
                        scoreMatrix[i0, i1, 0] = scoreInsertMax;
                    }

                    // layer 1: matches
                    scoreDiagMax = double.MinValue;
                    double matchCost = kernel.ComparisonScore(i0, i1);
                    for (int iLayer = 0; iLayer < 3; iLayer++)
                    {
                        scoreDiag = scoreMatrix[i0 - 1, i1 - 1, iLayer] + matchCost;
                        if (scoreDiag > scoreDiagMax)
                        {
                            scoreDiagMax = scoreDiag;
                            traceBackMatrix[i0, i1, 1] = iLayer;
                        }
                    }

                    if (bothFivePrimeEndsFree)
                    {
                        jumpScore = cumulativeScoreQuery[i0 - 1] + matchCost;
                        if (jumpScore > scoreDiagMax)
                        {
                            scoreDiagMax = jumpScore;
                            traceBackMatrix[i0, i1, 1] = -1;
                        }
                    }

                    scoreMatrix[i0, i1, 1] = scoreDiagMax;

                    // layer 2: deletions from sequence 0
                    scoreDeleteMax = double.MinValue;
                    for (int iLayer = 1; iLayer < 3; iLayer++)
                    {
                        scoreDelete = scoreMatrix[i0, i1 - 1, iLayer] + (iLayer == 2 ? kernel.ContinueDeletionScore(i0, i1) : kernel.OpenDeletionScore(i0, i1));
                        if (scoreDelete > scoreDeleteMax)
                        {
                            scoreDeleteMax = scoreDelete;
                            traceBackMatrix[i0, i1, 2] = iLayer;
                        }
                        scoreMatrix[i0, i1, 2] = scoreDeleteMax;
                    }
                }
            }

            return GetTracebackStart(qThreePrimeEndFree, rThreePrimeEndFree);
        }

        private int[] GetTracebackStart(bool qThreePrimeEndFree, bool rThreePrimeEndFree)
        {
            // find starting point and alignment score for traceback
            double maxScore = double.MinValue;
            int i0AtEnd = -1;
            int i1AtEnd = -1;
            int lastQIndex = dimQuery - 1;
            int lastRIndex = dimReference - 1;

            if (qThreePrimeEndFree && rThreePrimeEndFree)
            {
                double currentScore = 0;
                maxScore = currentScore;
                for (int i0 = 0; i0 < dimQuery; i0++)
                {
                    for (int i1 = 0; i1 < dimReference; i1++)
                    {
                        currentScore = scoreMatrix[i0, i1, 1] - cumulativeScoreQuery[i0]; 
                        if (currentScore > maxScore)
                        {
                            maxScore = currentScore;
                            i1AtEnd = i1;
                            i0AtEnd = i0;
                        }
                    }
                }
                AlignmentScore = maxScore + cumulativeScoreQuery[lastQIndex];
            }
            else if (qThreePrimeEndFree) // the traceback must start at the end of the reference sequence
            {
                double currentScore = 0;
                i1AtEnd = lastRIndex;
                maxScore = currentScore;
                for (int i0 = 0; i0 < dimQuery; i0++)
                {
                    currentScore = scoreMatrix[i0, lastRIndex, 1] - cumulativeScoreQuery[i0];
                    if (currentScore > maxScore)
                    {
                        maxScore = currentScore;
                        i0AtEnd = i0;
                    }
                }
                AlignmentScore = maxScore + cumulativeScoreQuery[lastQIndex];
            }
            else if (rThreePrimeEndFree)// traceback starts from end of query
            {
                double currentScore = 0;
                maxScore = currentScore;
                i0AtEnd = lastQIndex;
                for (int i1 = 0; i1 < dimReference; i1++)
                {
                    currentScore = scoreMatrix[lastQIndex, i1, 1];
                    if (currentScore > maxScore)
                    {
                        maxScore = currentScore;
                        i1AtEnd = i1;
                    }
                }
                AlignmentScore = maxScore;
            }
            else
            {
                i0AtEnd = lastQIndex;
                i1AtEnd = lastRIndex;
                maxScore = scoreMatrix[lastQIndex, lastRIndex, 1];
            }
            return new int[] { i0AtEnd, i1AtEnd };
        }

        /// <summary>
        /// Performs traceback for a dynamic programming-based pairwise alignment.
        /// </summary>
        public PolymerPair<T, U> TraceBack(int[] start, bool complete0 = false, bool complete1 = true)
        {
            // traceback either to the beginning, indicated by -1 as the next layer.

            PolymerPair<T, U> pair = new PolymerPair<T, U>(kernel.GetQuery(), kernel.GetReference());
            pair.AlignmentScore = this.AlignmentScore;

            if (start[0] < 0 || start[1] < 0)
            {
                pair.Index = new int[0][];
                pair.Score = new double[0];
                return pair;
            }

            int x = start[0];
            int y = start[1];
            int currentLayer = 1;
            if (start.Length > 2)
            {
                currentLayer = start[2];
            }

            Stack<int> indexXStack = new Stack<int>();
            Stack<int> indexYStack = new Stack<int>();
            Stack<double> scoreStack = new Stack<double>();

            int indexLength = 0;
            int nextLayer;
            bool finishing0 = false;
            bool finishing1 = false;
            while (currentLayer > -1)
            {
                nextLayer = traceBackMatrix[x, y, currentLayer];
                scoreStack.Push(scoreMatrix[x, y, currentLayer]);
                indexLength++;
                switch (currentLayer)
                {
                    case 0: // insertion
                        indexXStack.Push(x);
                        indexYStack.Push(-1);
                        x--;
                        break;
                    case 1: // diagonal move
                        indexXStack.Push(finishing1 ? -1 : x);
                        indexYStack.Push(finishing0 ? -1 : y);
                        if (x == 0) // reached top edge
                        {
                            finishing1 = true;
                            if (!complete1) // end the alignment here
                            {
                                nextLayer = -1;
                            }
                            y--;
                        }
                        else if (y == 0) // reached left edge
                        {
                            finishing0 = true;
                            if (!complete0) // end the alignment here
                            {
                                nextLayer = -1;
                            }
                            x--; 
                        }
                        else
                        {
                            x--;
                            y--;
                        }
                        break;
                    case 2: // deletion
                        indexXStack.Push(-1);
                        indexYStack.Push(y);
                        y--;
                        break;
                }
                currentLayer = nextLayer;
            }

            // declare the index array and fill it 
            pair.Index = new int[indexLength][];
            pair.Score = new double[indexLength];

            for (int i = 0; i < indexLength; i++)
            {
                pair.Index[i] = new int[2] { indexXStack.Pop(), indexYStack.Pop() };
                pair.Score[i] = scoreStack.Pop();
            }

            return pair;
        }
    }


    /// <summary>
    /// Provides the means to perform pairwise alignment.
    /// </summary>
    /// <typeparam name="T">The nucleotide type for the first polynucleotide.</typeparam>
    /// <typeparam name="U">The nucleotide type for the second polynucleotide.</typeparam>
    public class PairwiseAlignerSymmetric<T>
    {
        /// <summary>
        /// The total score for the alignment as reported.
        /// </summary>
        public double AlignmentScore { get; protected set; }

        private double[, ,] scoreMatrix;
        private int[, ,] traceBackMatrix;

        private double[][] priorScore;

        private int[] dimPhysical = new int[]{-1, -1};
        private int[] dim = new int[2];

        private AlignmentKernelSymmetric<T> kernel;

        /// <summary>
        /// Sets the default scores for the first column.
        /// </summary>
        /// <param name="value">The score to be accumulated over rows as the alignment scores that occur if the polynucleotide local alignment begins on a row other than the zeroth.</param>
        public void SetPriorScore(int sequenceNumber)
        {
            if (dim[sequenceNumber] < 1)
                return;

            for (int i = 0; i < dim[sequenceNumber]; i++)
            {
                priorScore[sequenceNumber][i] = kernel.PriorScore(sequenceNumber,i);
            }
        }

        /// <summary>
        /// Sets the default scores for the first column.
        /// </summary>
        /// <param name="value">The score to be accumulated over rows as the alignment scores that occur if the polynucleotide local alignment begins on a row other than the zeroth.</param>
        public void SetPriorScore(int iSeq, double score)
        {
            if (dim[iSeq] < 1)
                return;

            for (int i = 0; i < dim[iSeq]; i++)
            {
                priorScore[iSeq][i] = score;
            }
        }

        /// <summary>
        /// Sets the default scores for the first row.
        /// </summary>
        /// <param name="scores">A vector of alignment scores taken as the alignment scores that occur if the polynucleotide local alignment begins on a column other than the zeroth.</param>
        public void SetPriorScore(int iSeq, double[] scores)
        {
            if (priorScore[iSeq].Length != scores.Length)
            {
                throw new ArgumentException("Scores must have the same length as the existing prior score array.");
            }
            priorScore[iSeq] = scores;
        }

        /// <summary>
        /// Constructor.
        /// </summary>
        /// <param name="_pair">The polynucleotide pair to be aligned.</param>
        /// <param name="_comparator">A comparator of the appropriate types for use in the alignment.</param>
        public PairwiseAlignerSymmetric(AlignmentKernelSymmetric<T> kernel)
        {
            this.kernel = kernel;
            for (int iSeq = 0; iSeq < 2; iSeq++)
            {
                if (kernel.GetSeq(iSeq) == null)
                {
                    dim[iSeq] = dimPhysical[iSeq] = 0;
                }
                else
                {
                    dim[iSeq] = dimPhysical[iSeq] = kernel.GetSeq(iSeq).Seq.Length;
                }
            }

            // allocate matrices
            scoreMatrix = new double[dimPhysical[0], dimPhysical[1], 3];
            traceBackMatrix = new int[dimPhysical[0], dimPhysical[1], 3];
            priorScore = new double[2][];
            for (int iSeq = 0; iSeq < 2; iSeq++ )
            {
                priorScore[iSeq] = new double[dimPhysical[iSeq]];
                SetPriorScore(iSeq);
            }
        }

        /// <summary>
        /// Constructor.
        /// </summary>
        /// <param name="kernel">A comparator of the appropriate types for use in the alignment.</param>
        /// <param name="_dim0">The initial first dimension of the score and traceback matrices.</param>
        /// <param name="_dim1">The initial second dimension of the score and traceback matrices.</param>
        public PairwiseAlignerSymmetric(int dim0, int dim1, AlignmentKernelSymmetric<T> kernel)
        {
            this.kernel = kernel;
            this.dim[0] = dimPhysical[0] = dim0;
            this.dim[1] = dimPhysical[1] = dim1;

            // allocate matrices
            scoreMatrix = new double[dimPhysical[0], dimPhysical[1], 3];
            traceBackMatrix = new int[dimPhysical[0], dimPhysical[1], 3];
            priorScore = new double[2][];
            for (int iSeq = 0; iSeq < 2; iSeq++)
            {
                priorScore[iSeq] = new double[dimPhysical[iSeq]];
            }
        }

        public void SetSequence(Polymer<T> p, int sequenceNumber)
        {
            kernel.SetSeq(p, sequenceNumber);
            dim[sequenceNumber] = p.Seq.Length;
            if (dim[sequenceNumber] > dimPhysical[sequenceNumber])
            {
                dimPhysical[sequenceNumber] = dim[sequenceNumber];
                scoreMatrix = new double[dimPhysical[0], dimPhysical[sequenceNumber], 3];
                traceBackMatrix = new int[dimPhysical[0], dimPhysical[1], 3];
                priorScore[sequenceNumber] = new double[dimPhysical[sequenceNumber]];
            }
            SetPriorScore(sequenceNumber);
        }

        /// <summary>
        /// Fills the score matrix for pairwise alignment between the two Polynucleotides in the pair.
        /// </summary>

        /// <summary>
        /// Fills the score matrix for pairwise alignment between the two Polynucleotides in the pair.
        /// </summary>
        public int[] FillScoreMatrix()
        {
            dim[0] = kernel.GetSeq(0).Seq.Length;
            dim[1] = kernel.GetSeq(1).Seq.Length;

            int i0 = 0;
            int i1 = 0;

            double scoreDiag, scoreDiagMax;
            double scoreDelete, scoreDeleteMax;
            double scoreInsert, scoreInsertMax;

            // compute element [0,0]
            scoreMatrix[0, 0, 0] = priorScore[0][0];
            traceBackMatrix[0, 0, 0] = -1;
            scoreMatrix[0, 0, 1] = kernel.ComparisonScore(0, 0);
            traceBackMatrix[0, 0, 1] = -1;
            scoreMatrix[0, 0, 2] = priorScore[1][0];
            traceBackMatrix[0, 0, 2] = -1;

            double matchCost = double.MinValue;

            // compute zeroth row
            for (i1 = 1; i1 < dim[1]; i1++)
            {
                scoreMatrix[0, i1, 0] = double.MinValue;
                scoreMatrix[0, i1, 1] = kernel.ComparisonScore(0, i1);
                traceBackMatrix[0, i1, 1] = 2;
                scoreMatrix[0, i1, 2] = scoreMatrix[0,i1-1,2] + priorScore[1][i1];
                traceBackMatrix[0, i1, 2] = 2;
            }

            // compute zeroth column
            for (i0 = 1; i0 < dim[0]; i0++)
            {
                scoreMatrix[i0, 0, 0] = scoreMatrix[i0-1,0,0] + priorScore [0][i0];
                traceBackMatrix[i0, 0, 0] = 0;
                scoreMatrix[i0, 0, 1] = kernel.ComparisonScore(i0, 0);
                traceBackMatrix[i0, 0, 1] = 0;
                scoreMatrix[i0, 0, 2] = double.MinValue;
            }

            // compute entries up to last row and column
            for (i0 = 1; i0 < dim[0]-1; i0++)
            {
                for (i1 = 1; i1 < dim[1]-1; i1++)
                {
                    // layer 0: insertions into sequence 0
                    scoreInsertMax = double.MinValue;
                    for (int iLayer = 0; iLayer < 2; iLayer++)
                    {
                        scoreInsert = scoreMatrix[i0 - 1, i1, iLayer] + (iLayer == 0 ? kernel.ContinueGapScore(i0, i1) : kernel.OpenGapScore(i0, i1));
                        if (scoreInsert > scoreInsertMax)
                        {
                            scoreInsertMax = scoreInsert;
                            traceBackMatrix[i0, i1, 0] = iLayer;
                        }
                        scoreMatrix[i0, i1, 0] = scoreInsertMax;
                    }

                    // layer 1: matches
                    scoreDiagMax = double.MinValue;
                    matchCost = kernel.ComparisonScore(i0, i1);
                    for (int iLayer = 0; iLayer < 3; iLayer++)
                    {
                        scoreDiag = scoreMatrix[i0 - 1, i1 - 1, iLayer] + matchCost;
                        if (scoreDiag > scoreDiagMax)
                        {
                            scoreDiagMax = scoreDiag;
                            traceBackMatrix[i0, i1, 1] = iLayer;
                        }
                    }

                    scoreMatrix[i0, i1, 1] = scoreDiagMax;

                    // layer 2: deletions from sequence 0
                    scoreDeleteMax = double.MinValue;
                    for (int iLayer = 1; iLayer < 3; iLayer++)
                    {
                        scoreDelete = scoreMatrix[i0, i1 - 1, iLayer] + (iLayer == 2 ? kernel.ContinueGapScore(i0, i1) : kernel.OpenGapScore(i0, i1));
                        if (scoreDelete > scoreDeleteMax)
                        {
                            scoreDeleteMax = scoreDelete;
                            traceBackMatrix[i0, i1, 2] = iLayer;
                        }
                        scoreMatrix[i0, i1, 2] = scoreDeleteMax;
                    }
                }
            }

            // last row
            i0 = dim[0] - 1;
            for (i1 = 1; i1 < dim[1] - 1; i1++)
            {
                // layer 0: insertions into sequence 0
                scoreInsertMax = double.MinValue;
                for (int iLayer = 0; iLayer < 2; iLayer++)
                {
                    scoreInsert = scoreMatrix[i0 - 1, i1, iLayer] + (iLayer == 0 ? kernel.ContinueGapScore(i0, i1) : kernel.OpenGapScore(i0, i1));
                    if (scoreInsert > scoreInsertMax)
                    {
                        scoreInsertMax = scoreInsert;
                        traceBackMatrix[i0, i1, 0] = iLayer;
                    }
                    scoreMatrix[i0, i1, 0] = scoreInsertMax;
                }

                // layer 1: matches
                scoreDiagMax = double.MinValue;
                matchCost = kernel.ComparisonScore(i0, i1);
                for (int iLayer = 0; iLayer < 3; iLayer++)
                {
                    scoreDiag = scoreMatrix[i0 - 1, i1 - 1, iLayer] + matchCost;
                    if (scoreDiag > scoreDiagMax)
                    {
                        scoreDiagMax = scoreDiag;
                        traceBackMatrix[i0, i1, 1] = iLayer;
                    }
                }

                scoreMatrix[i0, i1, 1] = scoreDiagMax;

                // layer 2: deletions from sequence 0
                scoreDeleteMax = double.MinValue;
                for (int iLayer = 1; iLayer < 3; iLayer++)
                {
                    scoreDelete = scoreMatrix[i0, i1 - 1, iLayer] + priorScore[1][i1];
                    if (scoreDelete > scoreDeleteMax)
                    {
                        scoreDeleteMax = scoreDelete;
                        traceBackMatrix[i0, i1, 2] = iLayer;
                    }
                    scoreMatrix[i0, i1, 2] = scoreDeleteMax;
                }
            }

            // last column
            i1 = dim[1] - 1;
            for (i0 = 1; i0 < dim[0] - 1; i0++)
            {
                // layer 0: insertions into sequence 0
                scoreInsertMax = double.MinValue;
                for (int iLayer = 0; iLayer < 2; iLayer++)
                {
                    scoreInsert = scoreMatrix[i0 - 1, i1, iLayer] + priorScore[0][i0];
                    if (scoreInsert > scoreInsertMax)
                    {
                        scoreInsertMax = scoreInsert;
                        traceBackMatrix[i0, i1, 0] = iLayer;
                    }
                    scoreMatrix[i0, i1, 0] = scoreInsertMax;
                }

                // layer 1: matches
                scoreDiagMax = double.MinValue;
                matchCost = kernel.ComparisonScore(i0, i1);
                for (int iLayer = 0; iLayer < 3; iLayer++)
                {
                    scoreDiag = scoreMatrix[i0 - 1, i1 - 1, iLayer] + matchCost;
                    if (scoreDiag > scoreDiagMax)
                    {
                        scoreDiagMax = scoreDiag;
                        traceBackMatrix[i0, i1, 1] = iLayer;
                    }
                }

                scoreMatrix[i0, i1, 1] = scoreDiagMax;

                // layer 2: deletions from sequence 0
                scoreDeleteMax = double.MinValue;
                for (int iLayer = 1; iLayer < 3; iLayer++)
                {
                    scoreDelete = scoreMatrix[i0, i1 - 1, iLayer] + (iLayer == 2 ? kernel.ContinueGapScore(i0, i1) : kernel.OpenGapScore(i0, i1));
                    if (scoreDelete > scoreDeleteMax)
                    {
                        scoreDeleteMax = scoreDelete;
                        traceBackMatrix[i0, i1, 2] = iLayer;
                    }
                    scoreMatrix[i0, i1, 2] = scoreDeleteMax;
                }
            }

            // last entry
            i0 = dim[0] - 1;
            i1 = dim[1] - 1;
            // layer 0: insertions into sequence 0
            scoreInsertMax = double.MinValue;
            for (int iLayer = 0; iLayer < 2; iLayer++)
            {
                scoreInsert = scoreMatrix[i0 - 1, i1, iLayer] + priorScore[0][i0];
                if (scoreInsert > scoreInsertMax)
                {
                    scoreInsertMax = scoreInsert;
                    traceBackMatrix[i0, i1, 0] = iLayer;
                }
                scoreMatrix[i0, i1, 0] = scoreInsertMax;
            }

            // layer 1: matches
            scoreDiagMax = double.MinValue;
            matchCost = kernel.ComparisonScore(i0, i1);
            for (int iLayer = 0; iLayer < 3; iLayer++)
            {
                scoreDiag = scoreMatrix[i0 - 1, i1 - 1, iLayer] + matchCost;
                if (scoreDiag > scoreDiagMax)
                {
                    scoreDiagMax = scoreDiag;
                    traceBackMatrix[i0, i1, 1] = iLayer;
                }
            }

            scoreMatrix[i0, i1, 1] = scoreDiagMax;

            // layer 2: deletions from sequence 0
            scoreDeleteMax = double.MinValue;
            for (int iLayer = 1; iLayer < 3; iLayer++)
            {
                scoreDelete = scoreMatrix[i0, i1 - 1, iLayer] + priorScore[1][i1];
                if (scoreDelete > scoreDeleteMax)
                {
                    scoreDeleteMax = scoreDelete;
                    traceBackMatrix[i0, i1, 2] = iLayer;
                }
                scoreMatrix[i0, i1, 2] = scoreDeleteMax;
            }

            return GetTracebackStart();
        }

        private int[] GetTracebackStart()
        {
            // find starting point and alignment score for traceback

            int lastQIndex = dim[0] - 1;
            int lastRIndex = dim[1] - 1;
            int startLayer = 0;
            double maxScore = scoreMatrix[lastQIndex, lastRIndex, 0];
            for (int layer = 0; layer < 3; layer++ )
            {
                if (scoreMatrix[lastQIndex, lastRIndex, layer] > maxScore)
                {
                    startLayer = layer;
                    maxScore = scoreMatrix[lastQIndex, lastRIndex, layer];
                }
            }
            AlignmentScore = maxScore;
            return new int[] { lastQIndex, lastRIndex, startLayer };
        }

        /// <summary>
        /// Performs traceback for a dynamic programming-based pairwise alignment.
        /// </summary>
        public PolymerPair<T, T> TraceBack(int[] start)
        {
            // traceback either to the beginning, indicated by -1 as the next layer.

            PolymerPair<T, T> pair = new PolymerPair<T, T>(kernel.GetSeq(0), kernel.GetSeq(1));
            pair.AlignmentScore = this.AlignmentScore;

            if (start[0] < 0 || start[1] < 0)
            {
                pair.Index = new int[0][];
                pair.Score = new double[0];
                return pair;
            }

            int x = start[0];
            int y = start[1];
            int currentLayer = 1;
            if (start.Length > 2)
            {
                currentLayer = start[2];
            }

            Stack<int> indexXStack = new Stack<int>();
            Stack<int> indexYStack = new Stack<int>();
            Stack<double> scoreStack = new Stack<double>();

            int indexLength = 0;
            int nextLayer;
            while (currentLayer > -1)
            {
                nextLayer = traceBackMatrix[x, y, currentLayer];
                scoreStack.Push(scoreMatrix[x, y, currentLayer]);
                indexLength++;
                switch (currentLayer)
                {
                    case 0: // insertion
                        indexXStack.Push(x);
                        indexYStack.Push(-1);
                        x--;
                        break;
                    case 1: // diagonal move
                        indexXStack.Push(x);
                        indexYStack.Push(y);
                        if (x == 0) // reached top edge
                        {
                            y--;
                        }
                        else if (y == 0) // reached left edge
                        {
                            x--;
                        }
                        else
                        {
                            x--;
                            y--;
                        }
                        break;
                    case 2: // deletion
                        indexXStack.Push(-1);
                        indexYStack.Push(y);
                        y--;
                        break;
                }
                currentLayer = nextLayer;
            }

            // declare the index array and fill it 
            pair.Index = new int[indexLength][];
            pair.Score = new double[indexLength];

            for (int i = 0; i < indexLength; i++)
            {
                pair.Index[i] = new int[2] { indexXStack.Pop(), indexYStack.Pop() };
                pair.Score[i] = scoreStack.Pop();
            }

            return pair;
        }
    }
}
